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^ ! 1. Introduction 



ABSTRACT 



We define and put at the disposal of the community SOAP, Spot Oscillation And Planet, a software tool that simulates the effect of 
stellar spots and plages on radial velocimetry and photometry. This paper describes the tool release and provides instructions for its 
use. We present detailed tests with previous computations and real data to assess the code's performance and to validate its suitability. 
We characterize the variations of the radial velocity, line bisector, and photometric amplitude as a function of the main variables: 
projected stellar rotational velocity, filling factor of the spot, resolution of the spectrograph, linear limb-darkening coefficient, latitude 
of the spot, and inclination of the star. Finally, we model the spot distributions on the active stars HD166435, TW Hya and HD189733 
which reproduce the observations. We show that the software is remarkably fast, allowing several evolutions in its capabilities that 
could be performed to study the next challenges in the exoplanetary field connected with the stellar variability. 

Key words, methods: numerical - planetary systems - techniques: radial velocimetry, photometry - stellar activity - starspots 
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Nowadays, more than 700 exoplanet candidates are reported in 
the literature, most of them detected via the radial-velocimetry 
(RV) techniqu^H- However, this efficient method is indirect as 
well as the photometric transit or astrometry techniques. One of 
the problems with this is the fact that RV variations can in some 
cases be caused by other mechanisms that are not related to the 
presence of low-mass companions. Phenomena such as stellar 
pulsation, inhomogeneous convection, spots, or magnetic cycles 
can prevent us from finding planets, they might degrade the pa- 
rameters estimation, or give us false candidates, if they produce 
a stable periodic signal (e.g. Queloz et al. 2001; Bonfils et al. 
2007; Huelamo et al. 2008). An essential work is then needed 
to understand all phenomena caused by stellar variability, and 
to characterize their impact on RV and photometry to be able to 
distinguish between these and real planetary signatures. 

Particularly, dark spots and bright plages are present on the 
surface of dwarf stars from spectral types F to M, even in their 
low active phase (like the Sun, e.g. Meunier et al. 2010). Their 
appearance and disappearance on the stellar photosphere, com- 
bined with the stellar rotation, may lead to errors and uncertain- 
ties in the characterization of planets both in RV and photometry 
(e.g. Boisse et al. 2009). In a longer term, the change in the den- 
sity of spots along the stellar magnetic cycle can also induce 
variations in RV and spectroscopic indices (Santos et al. 2010; 
Dumusque et al. 201 1; Lovis et al. 201 1). 



* The tool is available at http://www.astro.up.pt/soap 

1 Consult the up-to-date catalogue of the website exoplanet. eu 



Before the detection of the first planet (Mayor & Queloz 
1995), it was already noticed that the rotation of spots with the 
stellar surface alters the line profiles, modifying the line cen- 
troids without changing the true RV of the star. It was also 
thought that line bisector variations are also induced by modifi- 
cation of the convection pattern. In 1997, Saar & Donahue devel- 
oped simple models from a medium-strength Fe I line at about 
6000A to examine the effect of starspots and convective inho- 
mogeneities. The authors derived a relation between the semi- 
amplitude of the RV, the v sin i of the star and the fraction of the 
stellar disk covered by the spot. Using models of Ca I 6439 A 
and Fe I 6430 A lines, Hatzes (1999, 2002) reproduced the RV 
and the bisector line deformation induced by spotted line profile 
and derived similar relations. Desort et al. (2007) computed the 
visible spectra of various spectral type stars and simulated the 
weight of the spot in each wavelength, applying a black-body 
law to the stellar spectrum as a function of the spot tempera- 
ture. To explain the variability observed in LkCa 19, Huerta et al. 
(2008) modeled the effect of starspots on the RV of young stars. 
To do so, they synthesized different spectra on a 100 A bandpass 
near 6300A with the spot depicted as a zone with a lower temper- 
ature than the stellar surface, with lower intensity and different 
spectral features. In a different way, Lanza et al. (2010) derived 
a model to compute the RV variation from the photometric light 
curve that they applied on CoRoT-7 and HD 189733 (Lanza et al. 
201 1). Aigrain et al. (201 1) presented a similar model. 

Stabilized spectrographs give high-precision RV measure- 
ments from the cross-correlation function (CCF) of the spectrum 
with a numerical mask (Baranne et al. 1996, Pepe et al. 2002). 
The CCF is then fitted with a Gaussian to derive the stellar RV. 
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The FWHM and contrast of the CCF were also recently found 
to be sensitive to short-term (Boisse et al. 2009; Queloz et al. 
2009) and long-term (Santos et al. 2010; Lovis et al. 2011) ac- 
tivity. 

This paper presents SOAP, a tool offered to the community 
that enables simulating spots and plages on rotating stars and 
computes their impact on RV and photometric measurements. 
It also computes the main spectroscopic diagnostic to monitor 
stellar activity: the bisector span (BIS), calculated as defined by 
Queloz et al. (2001), and delivers a model of the CCF that is 
equivalent to the mean line of the spectrum. In particular, we 
highlight a novel method that enables a fast computation of the 
spot position and geometry. 

The next section describes the software, giving the platform, 
the language, a short explanation of the principle of the code, 
a complete description of the input and output parameters, and 
finally, some details of the computing. We expose in a third sec- 
tion different tests to compare the results with other simulations 
and show the relation between the variations in RV, BIS and pho- 
tometry with the available parameters of the code. We then apply 
SOAP to the data of HD166435 (Queloz et al. 2001), TW Hya 
(Huelamo et al. 2008; Donati et al. 201 1) and HD189733 (Boisse 
et al. 2009). The last section discusses possible applications and 
future improvements that can be done with SOAP. 

2. Software description 

2.1. Platform and language 

The SOAP software is available for use at 
http://www.astro.up.pt/soap along with a brief manual 
and some example data. When using SOAP in publication, it is 
appropriate to cite this paper. 

SOAP is written with both C and python, for the core 
and high-level functions, respectively. An input parameter file 
driver . cf g has to be filled in by the user and uploaded 
on the webpage. After running the program, an output file, 
output . dat, and an output directory, CCFs, are returned. 

2.2. Principle 

SOAP first computes the non-spotted emission of the star in pho- 
tometry and in RV. To do that, the visible stellar disk is cen- 
tered on a grid of several tens of thousands of resolution ele- 
ments. The typical line of the emerging spectrum of the star is 
modeled by a Gaussian for each grid cell that is equivalent to 
the spectrum's CCF. Each Gaussian (in a given stellar position) 
is Doppler-shifted according to the projected rotational velocity 
and weighted by a linear limb-darkening law. 

SOAP calculates the position of the surface inhomogeneities 
defined by their latitudes, longitudes, and sizes. The cells inside 
the spots are modeled by the same Gaussian as for the stellar 
disk, but are weighted by their brightness (0 for a dark spot, 
higher than 1 for a bright plage). The code then removes (for 
dark spots) or adds (for bright spots or plages) the CCF and flux 
of the inhomogeneities to those of the non-spotted star. Finally, 
SOAP delivers the integrated spectral line, the flux, the RV, and 
the BIS as a function of the stellar rotational phase. 

2.3. Input parameters 

The input parameters are assigned by the user in the file 
driver, cfg. They are displayed in Table Q] and described in 
the following. 




Fig. 1. In three dimensions, the spot is first located in front of 
the observer (the x-axis is aligned with the line of sight), nrho 
number of points are equally distributed along the circumference 
of the spot (small black dots). 



Resolution of the simulation: The parameter grid determines 
the number of cells to resolve the stellar disk on a squared 
box of gridxgrid resolution elements. The position of the 
spot at each phase is settled thanks to the determination of 
the location of its circumference (cf. Sect. 3. 2. 2). The spot 
circumference is determined with nrho number of resolution 
elements (cf. Fig [TJ. We find that there are no significant 
changes in the results for values beyond grid-300 and 
nrho-2Q, and conclude that this option combines fast 
computation and sufficient resolution. 

Stellar parameters: The user can change the following stellar 
parameters: linear limb-darkening coefficient limb, stellar 
radius radius expressed in Sun radius, rotation period Prot 
in day, stellar inclination angle with the line of sight /, mean 
star's velocity gamma in kms~', and the initial phase for the 
computation of the simulation, psi. We note that SOAP fixes 
the reference of 0° in longitude as the position of the first 
spot in front of the line of sight at phase equals to 0. 

Spot parameters: The code allows including up to four spots, 
and choosing their longitude long and latitude lat in degrees. 
Their radius is defined by size and expressed as a linear 
projected value of the stellar radius. We remark that for 
a filling factor of less than 10%, a radius expressed as a 
linear projected value or as an arc value along the surface is 
completely identical. The brightness of the spot, defined by 
bright, corresponds to a weight to the stellar flux, and could 
take all positive values. Set to 0< bright < 1, it will model 
a darker spot than the stellar photosphere and, with a value 
higher than 1, it will simulate plages with a higher brightness 
than the star. This weight is applied to the CCF that modeled 
the emitted spectrum. We note that the same Gaussian is 
used for the photosphere and the inhomogeneities. 

CCF parameters: The parameters of the Gaussian can also be 
chosen to account for the resolution and accuracy of the 
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Fig. 2. Gaussian that modeled the CCF. The contrast of the CCF 
corresponds to the amplitude of the Gaussian. The window in RV 
to compute the CCF is the x-axis, here the window parameter is 
equal to 20 kms -1 . The red points illustrate the RV where the 
CCF values are calculated. The span between two points is the 
step parameter, here equal to 0.5 kms -1 . 

Table 1. Input parameters 
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other parameters, window and the step in kms" 1 , give the 
characteristics to calculate the CCF function. For low -rotator 
stars and high-resolution spectrographs, typical values are 
window=20kms~ l and sfep=0.5kms -1 (see Fig. [2]). Higher 
values should be used for window for young and/or fast- 
rotators stars. 
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-90 to +90 


SOAP simulates up to 4 spots 


size [RJ 




spot size in stellar radius 


bright 


Oto >1 


dark spot=0 ; plages>l 


profO 


0.4 




sigmaO [kms -1 ] 


2 to 5 


spectrograph resolution 


window [kms -1 ] 


20 


window calculating CCF function 


step [kms -1 ] 


0.1 





Fig. 3. Figure generated by SOAP (optional) showing the pho- 
tometry, the RV and the BIS (from top to bottom panel) as a 
function of the stellar rotational phase. As a default parameter, 
the first spot is in front of the line of sight at phase 0. 



be directly plotted by SOAP as a function of the stellar rotational 
phase if screen output is selected (Fig. [3). 

The phase step can be chosen in two ways. It can be a con- 
stant step of a fraction of the phase, phstep between 0. to 1. 
If this last is set equal to 0, the code will read phases in a file 
determined as an input parameter phJn. 

The computed CCFs can also be stored in the directory 
named CCFs if the ccfjout is set to 1 (no output = 0). All input 
parameters are archived in the header of the FITS files gener- 
ated. Each contains a one-dimension FITS file of the flux as a 
function of the RV on a window and with a step determined by 
the input parameters of the CCF. 



spectrograph, as well as the way the observed CCF was 
computed to be as close as possible to real data. The code 
asks for the width in kms -1 , sigmaO, of the CCF of a non- 
rotating star (or with rotational velocity too low to be re- 
solved) with the same spectral type. We recall that it is re- 
lated to the FWHM by cr = E^MM, in kms" 1 . The two 



2.5. Detailed numerical computing 
CCF and flux of the non-spotted star 

The star is computed as a disk with a radius normalized to 1. 
This disk is centered on a grid of grid x grid cells on a y-z plane. 
The cells take their y and z values between - 1 and 1 . Each grid 
cell (p v ,p z ) contains a flux value and a CCF. The flux value de- 
pends on the linear limb-darkening law. The CCF is modeled by 
a Gaussian with width and amplitude given by the input parame- 
ters, Doppler-shifted according to the projected rotational veloc- 
ity and weighted by the linear limb-darkening law. The projected 
rotational velocity is calculated from the rotational period, stel- 
lar radius and inclination of the spin axis with the line of sight, 
given in the input parameters. 

The code first sums the contribution of all grid cells of the 
non-spotted star to generate reference CCF and flux value. 



2.4. Output parameters 

The format of the output parameters is determined in the file 
driver . cfg. SOAP gives the values of the flux, RV and, BIS, as 
a function of the phase step and writes the result in the output file 
file_out. By default filejout is named 'output.dat' . It gives a file 
with four columns: phase, flux, RV, and BIS. These values can 



Position and contribution of the inhomogeneities ("spots") 

The spot is defined in the simulation as a disk on the stellar 
surface. When the spot rotates with the stellar photosphere, 
the shape of the spot seen along the line of sight changes 
and becomes ellipsoidal. The code is able to compute this so 
fast (see Sect. 13. U because it does not calculate the projected 
geometry of the spot. 
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The grid is on the y-z plane and the x-axis is along the line 
of sight. If the star is seen edge-on, it rotates around the z-axis. 
We remember that the stellar angle in the plane of the sky has 
no impact on the result and is completely omitted. The grid 
cells take values between -1 and 1. These maximum values 
correspond to the stellar radius. 

For each phase, the code first determines the area of the grid 
where the inhomogeneities are located. To do that, the spot is 
initialized in front of the line of sight, on the equator of a star 
seen edge-on. From there, the spot has a well known geometry 
of a circle centered on the x-axis and a radius size (in stellar 
radius) given in the input parameters. 

For each nrho, the number of points on the circumference of 
the spot, the code calculates the position according to the stellar 
inclination, spot latitude and longitude (see step 1 in Fig. @). 
For that a three-rotations matrix is applied to the nrho points: a 
rotation of the angle of the stellar inclination i about the y-axis 
Ry((f), a rotation of the spot longitude long about the z-axis R z {9), 
and a rotation of the spot latitude lat about the y-axis R y {if/)'. 



R y (<p)xR z (e)xR y ^) = 



with 



a — sin 9 sin if b 
sin 9 cos if/ cos 8 sin 8 sin if/ 
c sin if sin 9 d 



(1) 



a — cos ip cos 9 cos if/ - sin if/ sin if 
b = cos if cos 9 cos if/ + sin f cos if/ 
c = - sin if cos 8 cos if/ - cos if sin if/ 
d = cos if cos if/ - sin if cos 9 sin if/. 

We note that the conventions give tp -90 — i (in degree) and 
if/ = -lat and at that point, the xyz values are given as a fraction 
of the stellar radius. 

The xyz positions of the nrho points are then rotated accord- 
ing to the rotational phase, with an angle <p-360xph, ph is a 
number between and 1 (see Sectior l2~4l i. The rotation is per- 
formed along the unit vector aligned with the stellar inclination 
h=(cos(/), 0, sin(/)). The rotation matrix R u ((f>) is 



with 



mix + cos ffUxUy - /?u z cni x u z + /3\x y 
ffu y u x + /3u z mi y + cos <p cniyUz - /3u x 
oru z u x - ygUy mi z u y + /3u x au\ + cos <p 



(2) 



a - ( I - cos <p) 
j3 = sin <p 

From these nrho xyz positions, the code defines the visibility 
of the spot. If it is visible, it defines a smaller area where the 
spot is according to the minimum and maximum values taken by 
y and z (see step 2 in Fig. 0J. When there are both visible and 
invisible points, i.e. the spot is on the edge, the minima/maxima 
are over-/underestimated if the spot is on one of the axis (y or z). 
This is solved following the idea that if the spot is on the bottom 
z-axis (z < 0), then the minimum z value is equal to - 1 . 

Then, the maximal/minimal y-z values are transposed onto 
the grid (p y ,p z ), which depends on the sampling given by the 
grid parameter (see step 3 in Fig. |4j». We note that the code never 
memorizes the position of the spot in the grid coordinate (p y , p z ). 

Then the code calculates the impact of the spot on the stellar 
photometry and velocity. The code scans the smaller grid and 



checks whether each center of the cell is located inside or outside 
the spot. Because we do not know the projected geometry of the 
spot in its actual position, we perform an inverse rotation of the 
matrix given in EqQ] to replace the grid point where it would 
be in the initial configuration, an equatorial spot on an edge- 
on star, where the spot has the well-known geometry of a circle 
centered on the x-axis (see step 4 in Fig. @J. If the grid-point 
(p y ,p z ) is inside the spot, a Gaussian is attributed to this point 
with a velocity according to the projected rotational velocity, and 
weighted by the linear limb-darkening law and the intensity of 
the spot defined by its brightness (1 - bright). For this last point, 
we note that the computed CCF for the spots is to represent what 
is removed (for dark spots) or added (for bright spots) to the 
unspotted stellar CCF. The code then generates a summed CCF 
and a flux value for the spots at each wanted rotational phase. 
For dark spots, the CCF and flux values correspond to the stellar 
intensity that is hidden by the spots. 

Final parameters 

The code then subtracts for each phase the spots CCF from the 
non-spotted stellar CCF, and the spots flux value from the non- 
spotted stellar flux. SOAP returns a CCF and a flux value for 
each queried phase. 

The flux value corresponds to the photometry, which is nor- 
malized to the maximum value along the rotational phase. 

Each CCF is normalized to the maximum value of each CCF. 
Each CCF is fitted by a Gaussian to return RV values for each 
phase. The BIS values are calculated as described in Queloz et 
al. (2001). 



3. Tests and comparison with data 

3.1. Performance 

The computation of the code is sufficiently fast to allow many 
simulations and adjusting the data. For a MacBook Pro with a 
2.33 GHz Intel Core Duo, calculating 10,000 phases of a star 
with one always visible spot takes less than 40 s. To simulate 
10,000 phases of a star with four spots, the calculation and the 
output files are ready in less than 90 s. We emphasize that in 
the simulations reported below (Sect.4.3), a sufficient number 
of 100 phases per rotational period were computed. The code is 
therefore perfectly suitable for solving inverse problems . 

3.2. Looking at different parameters 

In Boisse et al. (2011), we used SOAP to derive the periodicity 
properties of the RV variations induced by spots, as well as to 
derive the relations between the RV, the photometry and the CCF 
parameters (BIS, FWHM and contrast). For the simulations, we 
chose a GOV star with a radius of 1.1 R and a vsin/=5.7kms _1 , 
a linear limb darkening of 0.6, and a spectrograph resolution of 
110 000. The spot was considered as a dark surface without any 
emission of light, and a size of 1% of the visible stellar surface. 
We refer to this paper for the computation of the RV variations 
induced by one or two spots with different stellar inclinations 
and spot latitudes and longitudes (cf. their Figures 1 and 11), 
and also for the corresponding relations between the RV and the 
bisector span (cf. their Figures 6 and 13). These authors also 
showed the relations between the FWHM and the photometry for 
one and two spots (their Figures 8 and 14). For one specific case, 
the relation between RV, BIS, FWHM, contrast, and photometry 
is shown in their Figures 1 and 7. On the webpage, one can find 
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Fig. 4. From left to right: 1) According to the stellar inclination and to the spot's latitude and longitude, the xyz position of the nrho 
points on the spot's circumference are calculated in 3-D. 2) The maximum and minimum y and z determine a smaller area where the 
spot is. 3) These \yz) positions (in 3D) are projected on the y-z plane, which gives (p y ,p : ) positions. 4) From these values, a smaller 
grid is defined. Each cell of this grid is scanned. A reverse to that performed in step 1 rotation is performed to establish if the cell is 
inside the spot or not. 
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Fig. 5. Semi-amplitude of the RV (squares) and BIS (circles) 
variations as a function of v sin i. The empty points are for a res- 
olution power of the spectrograph of R= 75,000 and the filled 
ones are for R= 40, 000. All other parameters are fixed (cf. text). 
The lines are the best chi-square fits from equations [3]and 2] 



BIS (circles) as a function of the v sin i. The semi-amplitude is 
defined in this paper as half of the maximum deviation, hence 
^yS, and later f° r tne photometry. We considered ob- 
servations made with a spectrograph with a resolution power of 
about R= 75,000 (empty points) and R= 40,000 (filled ones). 
We recall that the resolution power is coded in SOAP as sigmao, 
the width of the CCF of a non-rotating star or with a rotational 
velocity too low to be resolved. The sigmao values of GO star 
are derived from the Eq.B.2 and B.3 of Boisse et al. (2010) for 
a spectrograph resolution of 75,000 and 40,000. In contrast with 
previous results, we observe that the RV depends not only on the 
v sin i but also slightly on the resolution of the spectrograph. As 
shown by Desort et al. (2007), the BIS variations also depend 
on v sin i and on the resolution of the spectrograph. This depen- 
dence on the instrumental resolution is explained because the RV 
and the BIS are derived from the fit of a Gaussian on the CCF 
(cf. Section I2.51 l and the deformation of the CCF induced by a 
spot of a given size depending on the instrumental resolution. 

We derived the relation between the semi-amplitude of the 
variation of the RV, BIS, and filling factor of the spot f r , plotted 
in FigM 

We obtained the following laws: 



ARV 



ABIS 



20.67 f r ( Jv sin i 2 + sigma^ - sigmao) 9& 



9.6 



V sin i 2 + sigma 2 - sigmao) 



1.65 



(3) 



(4) 



computed values (output files) for specific cases so that the user 
can check the proper calculation of the code. 

In Fig. [3j we show the flux, RV, and BIS as a function of 
the stellar rotational phase for one equatorial spot of 1 % on the 
surface of an edge-on star with equivalent parameters that are 
given in the previous paragraph. 

To check the reliability of the code with previous simula- 
tions, we sought to reproduce the laws derived first by Saar & 
Donahue (1997), and by Hatzes (1999, 2002) and Desort et al. 
(2007) (cf. introduction). We simulated a spotted star as close as 
possible to the star of the Saar & Donahue (1997) simulation. 
We generated a 1 % dark spot on the equator of a G2 V edge-on 
star. In the following analysis of this section, all parameters are 
fixed on these values, except those studied in each paragraph. 
We plot in Fig.|5]the semi-amplitude of the RV (squares) and the 



These relations give an equivalent range of the amplitudes to 
those derived by Desort et al. (2007). The differences emerge be- 
cause we did not simulate a temperature for the spot. We remark 
that Desort et al.'s results are close to those of Saar & Donahue 
(1997) and Hatzes (1999, 2002) at low v sin; but then they de- 
part strongly. As Desort et al. (2007) noted, this is because Saar 
& Donahue (1997) and Hatzes (1999, 2002) modeled the effect 
in only one single line that has a varying sensibility to the impact 
of spots. Hatzes (1999) also noticed that. As already mentioned 
by Desort et al. (2007), these estimates, derived by us in Eqs. [3] 
and|4]orby the previously quoted authors, are indicative and we 
discourage quantitative conclusions derived from the blind use 
of one of these equations. 

In Fig. [5] we observe that when v sin i is lower than a cer- 
tain value, depending on the resolution of the spectrograph, the 
bisector does not vary, and thus cannot be used as a diagnostic 
of stellar activity. This phenomenon was already pointed out by 
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tometry (triangle) as a function of the linear limb-darkening co- 
efficient. All other parameters are fixed. 



Santos et al. (2003) and was already observed for active stars 
with a very low rotational velocity (e.g. Bonfils et al. 2007). 

We also aimed to characterize the impact of the linear limb- 
darkening coefficient. Following Claret (2004), this coefficient 
varies roughly from 0.5 to 0.9 (in V band) for dwarf stars with 
an effective temperature between 6500 K and 3800 K. Fig. [7] il- 
lustrates that for the same spot, a higher linear limb-darkening 
coefficient induces an increased impact in RV, BIS, and photom- 







Fig. 8. RV, BIS, and photometry as a function of the stellar rota- 
tional phase. The different colors illustrate the change in latitude 
of the spot. The maximum amplitude variation of the three pa- 
rameters are given for an equatorial spot. No variation is detected 
when the spot is on the pole. We note that beyond a certain value 
of lat (» 50°), the shape of the BIS variation changes. 



etry. This is mainly explained by a total lower brightness of the 
spotless star. 

Figs. l8|[T0l show the variations in the amplitude of the RV, 
BIS, and photometry as a function of stellar inclination and spot 
latitude angles. We remark that the variation laws are close to 
a cos 2 , as drawn with dashed lines in the figures, as expected 
because of the effects of projection. 

We emphasize that in our simple model, most of the variables 
are degenerate (like the size and the brightness). Nevertheless, 
some parameters can be constrained by observations that allow 
to model real data. 



3.3. Simulations of published spotted stars 

Several publications reported stellar RV variations that are 
caused by magnetic activity and not to a gravitationally bound 
companion. These cases could be reproduced by SOAP, which 
sheds light onto the parameters of the spotted stars. 



HD1 66435 

Queloz et al. (2001) published RV variations with a period of 
3.7987 days of HD166435 observed during two years. The star 
was found to be photometrically and magnetically variable with 
the same period. This and the long-term phase instability of the 
RV signal, are signatures of the evolution of a spot in a sta- 
ble active region. We simulated the HD 166435 variability with 
one spot on the stellar surface. This GOV star was observed 
with ELODIE, a spectrograph with a resolution power of about 
42,000. Queloz et al. (1998) published for this spectrograph the 
relation between the stellar (B - V) and sigmao, the width of the 
CCF of a non-rotating star or with a rotational velocity too low 
to be resolved. We used their equation (2) and the HD 166435 
(B- V)=0.633 value to derive a sigmaO equal to 4.62 kms -1 . We 
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Fig. 9. Semi-amplitude of RV (square), BIS (circle), and pho- 
tometry (triangle) as a function of the spot latitude. All the other 
parameters are fixed. The dashed lines correspond to a cos(/) 2 fit. 
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Fig. 10. Semi-amplitude of RV (square), BIS (circle), and pho- 
tometry (green triangle) as a function of the inclination of the 
star with the line of sight i. Two different v sin i at 3 and 7 kms" 1 
(filled and empty markers) are drawn. The photometry does not 
depend on vsim. All other parameters are fixed (cf. text). The 
dashed lines correspond to a cos(/) 2 fit. 



chose a linear limb-darkening coefficient of 0.6. Queloz et al. 
(2001) derivedthe stellar v sin /=7.6±0.5kms _I from the width of 
the CCF. With a radius of 0.99R o (Takeda 2007), we can then de- 
rive a stellar inclination of 35° with the line of sight. The bright- 
ness and the size parameters are almost completely degenerate 
for small spots. Therefore, we decided to fix bright to 0. and to 
adjust the size of the spot. Because only smooth variations of the 
RV, BIS, and flux (without plateau) are observed, this suggests 
that the spot is always visible to the observer. Using the sim- 
ulations, smooth variations of the parameters are observed for a 
spot with a latitude greater than 45 ° at a stellar inclination of 35 . 



~ 0.99 



u 0.98 




0.4 0.6 0.8 

Phase 



Fig. 11. SOAP simulation of a star with the same characteristics 
as HD166435 observed with ELODIE with a spot of 2.02% of 
the stellar surface to reproduce the RV, BIS, and flux variations 
observed by Queloz et al. (2001). 



Finally, we varied the latitude and size of the spot to be as close 
as possible to the observations. To obtain a good ratio between 
the RV and BIS semi-amplitude, we need to have a spot latitude 
close to 68°. The nearest semi-amplitude values to the observa- 
tions are then obtained for a spot size of 2.02% of the stellar 
disk. The related normalized flux varies by 2.9%, which corre- 
sponds for a relative magnitude of 6.85, to 0.033mag, which is 
in the lower range of the observed 0.035 to 0.05 mag variation. If 
we fit the RV variations with a Keplerian model, an eccentricity 
of 0.22 is found, equivalent to the 0.2 value reported by Queloz 
et al. (2001). We then agree with the conclusion of Queloz et 
al. (2001) that a simple one-spot model can well reproduce the 
observed variability. In Fig. [TT] we show the result of our simu- 
lation with the input parameters given in Table [2] We tested if a 
two-spot scenario at opposite longitudes could also explain the 
data. The period would then be equal to 7.5974d. With a radius 
of O.99R , the v sin i should be overestimated. We chose to put 
the star equator-on, which corresponds to a vsin/ of 6.6 kms 1 . 
We fixed the brightness of the spot to 0. We ran the simulations 
and found that no solution leads to the measured values of the 
semi-amplitude of the RV and the BIS (and to a good ratio be- 
tween the RV and BIS semi-amplitude variations). We conclude 
that the two-spot scenario is not appropriate. 



TW Hya 

The classical T Tauri star, TW Hya, was first announced to host 
a planet (Setiawan et al. 2008) before infrared RV measurements 
showed that a cold spot covering ~7% of the stellar surface and 
located at a latitude of 54° better explained the observations 
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Table 2. Input parameters for the simulation of HD 166435 



Parameter 


Value 


Ket. 


grid 


300 




nrho 


20 




limb 


0.6 


GOV star 


radius [R Q ] 


0.99 


Takeda 2007 


Prot [day] 


3.7987 


Queloz etal. 2001 


i[°] 


35 


calculated (see text) 


n -It 

gamma [kms J 


0. 




psi 


0. 




lat [°] 


68 




bright 


0. 


fixed 


size [R*] 


0.142 




sigmaO [kms -1 ] 


4.62 


ELODIE resolution 


window [kms -1 ] 


20 




step [kms -1 ] 


0.1 





Table 3. Input parameters for the simulation of TWHya 



Parameter 


Value 


T > p 

Ket. 


grid 


300 




nrho 


20 




limb 


0.6 




radius [R G ] 


1.1 


Donati et al. 2011 


Prot [day] 


3.56 


Setiawan et al. 2008 


i[°] 


15 


Donati et al. 2011 


gamma [kms J 


0. 




psi 


0. 




lat [°] 


80 




bright 


0. 


fixed 


size [R*] 


0.14 




sigmaO [kms -1 ] 


4.6 


ESPaDOnS resolution 


window [kms -1 ] 


20 




step [kms -1 ] 


0.1 





(Huelamo et al. 2008). Indeed, the infrared RV variations pre- 
sented a lower amplitude than the optical ones, in agreement 
with the spot scenario. The contrast between the spot and the 
photosphere is weaker in the near-infrared than in optical, in- 
ducing weaker RV variations. The spot model was derived with 
SOAP using the optical observations. All parameters needed for 
the simulations are given in the Section 2 of Huelamo et al. 
(2008). These authors reproduced an observed semi-amplitude 
of about 150ms -1 . The model predicts a photometric variation of 
4%, much lower than the observed 20%, which prompted the au- 
thors, among other indices, to conclude that there is a dark spot 
combined with a bright one due to accretion. Recently, Donati 
et al. (201 1) observed TW Hya with the ESPaDOnS spectropo- 
larimeter at the CFHT. These authors confirmed that the RV fluc- 
tuations should come from a high-latitude photospheric cool spot 
overlapping with the main magnetic pole where most of the ac- 
cretion is concentrated. The authors observed peak-to-peak RV 
variations of a;50ms -1 that they assumed to originate from a spot 
at high-latitude with a size of 2% of the stellar surface. We used 
SOAP to simulate their data with the input parameters given in 
Table [3] As observed by Setiawan et al. (2008) and Huelamo et 
al. (2008), the computed BIS variations are too small to be de- 
tected by instruments with a resolution power lower than 70,000. 
The computed RV variations agree with the Donati et al. (201 1) 
observations. But the simulated photometric variability is even 
lower than in the simulation by Huelamo et al. (2008) with less 
than 1%. While the spot may have evolved between the different 
observations, we can rather conclude that a simple spot model is 
not sufficient to explain all variables, and a combination of a dark 
spot and bright accretion probably explains the observables of 
this young star. We note that a hot accretion spot should be seen 
as a continuum locally in the stellar photosphere and may also 
induce RV drift. However, because hot spots are usually smaller 
than cold ones, we expect their impact to be significantly weaker. 

HD1 89733 

HD 189733 is an active star that hosts a transiting hot-Jupiter. 
Boisse et al. (2009) reported ~2-months observations performed 
with the SOPHIE spectrograph simultaneously to photometric 
MOST data to monitor the stellar activity. We used their orbital 
solution to remove the effect of the hot Jupiter in the RV data and 
modeled the residuals. The simulation was performed with one 
dark spot and the input parameters reported in Table We chose 
a linear limb-darkening value of 0.8 according to Claret (2004) 



Table 4. Input parameters for the simulation of HD 189733 



Parameter 


Value 


Ref. 


grid 


300 




nrho 


20 




limb 


0.8 


K2V star 


radius [R Q ] 


0.766 


Triaud et al. 2009 


Prot [day] 


11.95 


Henry & Winn 2008 


i n 


85.7 


Pont et al. 2007 


gamma [kms -1 ] 


0. 




psi 


0. 




lat [°] 


70 




bright 


0. 


fixed 


size [R*] 


0.14 




sigmaO [kms -1 ] 


2.68 


SOPHIE resolution 


window [kms -1 ] 


20 




step [kms -1 ] 


0.1 





for [M/H]=0, T eff =5000K and logg=4.5. In Fig.7 of Boisse et 
al. (2009), we observed a plateau in the variations of the pa- 
rameters as a function of the rotational phase. This plateau is 
expected when the main spot is out of sight. We obtain a semi- 
amplitude close to the observed variations with a spot at 70° lat- 
itude and size of about 2% of the stellar disk, ARV/2=14ms -1 , 
ABIS/2=6ms -1 and AFlux of 1%. The derived size of the spot 
agrees with the value estimated by HST observations of the star 
by Pont et al. (2007), thanks to occultations of spots by the tran- 
siting planet. Even if we can expect several spots in the photo- 
sphere of this star (e.g. Pont et al. 2007; Lanza et al. 2011), we 
observed that a one-spot model can well explain the observed 
RV BIS, and flux variations (Fig. P. 

4. Perspectives 

At this stage, and as presented in this paper, SOAP is at a quite 
basic model level. Because it is fast and clearly coded, improve- 
ments and developments can be easily implemented. Dumusque 
et al. (2011) already used this tool to simulate the evolution of 
the number of spots along the magnetic cycle of the Sun, and in- 
ferred the corresponding RV amplitude. These authors deduced 
a method to average the effect of the stellar activity to search for 
low-mass planets. 

Additional improvements will be made on the tool. SOAP 
simulates the photometric effect due to inhomogeneities on the 
stellar surface. However, it does not consider the effect that they 
also block the convection pattern. Moreover, spots and plages 
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Phase 



Fig. 12. In blue lines, SOAP simulation of a star with the same 
characteristics as HD 189733 observed with SOPHIE with a spot 
of 2% of the stellar surface to reproduce the RV, BIS, and flux 
variations observed by Boisse et al. (2009) (black data points). 



are not only areas darker or brighter, but have a different tem- 
perature than the stellar photosphere. They then present differ- 
ent spectra than the stellar photosphere, which could be seen as 
different CCF parameters. This wavelength dependency must be 
taken into account because, if most of the accurate spectrographs 
are now in the optical, soon near-infrared instruments will be 
available, a wavelength domain in which the amplitude of RV 
variation due to activity is supposed to be smaller (e.g. Huelamo 
et al. 2008; Prato et al. 2008; Figueira et al. 2010; Crockett et al. 
2011;Mahmud etal. 2011). 

Stellar magnetic fields might be monitored by po- 
larimetry of the stellar lines induced by Zeeman ef- 
fect (Donati & Landstreet, 2009) with instruments such 
as ESPADONS @CFHT, Narval@TBL, HARPSPol@ESO, or 
soon SPIRou@CFHT also dedicated to the search of planets 
around low-mass stars. A first attempt to look for relations be- 
tween RV variations and polarized light is promising (Delfosse 
et al., in prep.). SOAP could then be developed to simulate the 
polarimetric emission of the active star. 

When the planets are transiting in front of their parent stars, 
their radii are measured. And when they are also characterized 
in RV, which yields the mass, it leads to a measure of the mean 
density, which is the first step in the internal characterization of 
exoplanets. But the stellar inhomogeneities may have impacts 
on the radius determination derived from transit detections. 
One of the impacts is that it hampers a correct determination 
of the radius. That will be crucial for determining atmospheric 
components that are derived from the radius variations as a 
function of wavelength (e.g. Czesla et al. 2009). 

Dedicated observations of active stars are therefore needed, 
if possible simultaneously, in RV, photometry, and polarimetry 



to characterize the different level of stellar variability, to under- 
stand the correlation between the different proxies, and to con- 
strain the simulations. The release of thousands of light curves 
from the space missions CoRoT and Kepler dedicated to astero- 
seismology and the search for planets via transit will help us to 
understand statistically the photometric variability of these stars. 

This tool is useful for computing and testing the effect of 
stellar activity on RV and photometry. It will help to understand 
the challenges related to the knowledge of stellar activity for the 
next decade: detect rocky planets in the habitable zone of their 
stars (from G to M dwarfs), understand the activity in the low- 
mass end of M dwarf (on which future projects such as SPIRou 
or CARMENES will focus), limitation to the sum of different 
transit observations to characterize the atmospheric components 
(from the ground or with Spitzer and JWST), and search for 
planets around young stars. These can be simulated with SOAP 
to search for indices and corrections of the effect of activity. 
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